Dynamical evolution of unstable self-gravitating scalar solitons 
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Recently, static and spherically symmetric configurations of globally regular self-gravitating scalar 
solitons were found. These configurations are unstable with respect to radial linear perturbations. 
In this paper we study the dynamical evolution of such configurations and show that, depending on 
the sign of the initial perturbation, the solitons either collapse to a Schwarzschild black hole or else 
"explode" into an outward moving domain wall. 
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I. INTRODUCTION 

Recently, a new family of scalar-hairy black holes (BH) 
and their corresponding solitons (scalarons) were found 
within an Einstein-Higgs theory with a non-positive semi- 
definite scalar field potential V{(j)) Q. This kind of po- 
tential violates the weak-energy condition (WEC) and 
therefore invalidates the applicability of the no-scalar- 
hair theorems These configurations are interest- 
ing in several respects. On one hand, they constitute an 
example that obstructs the extension of no-hair theorems 
to potentials of this type. On the other hand, they can 
be useful for testing some of the predictions of the recent 
isolated-horizons formalism In fact, these configura- 
tions can be shown to be unstable with respect to radial- 
linear perturbations, and therefore they can be seen as 
bound states of non- hairy black holes and scalarons (cf . 
in the context of colored BH) . The simple perturbation 
analysis, however, does not provide any definite answer 
about the final fate of these configurations. Neverthe- 
less, an heuristic analysis based on energetic arguments 
does provide some clues about their fate. Presumably, 
the plain Schwarzschild BH constitutes the lower energy- 
mass bound (the "ground state" ) of possible BH configu- 
rations with fixed boundary conditions, which correspond 
to fixed horizon area Ah and asymptotic flatness. There- 
fore, among all BH configurations within the theory, the 
Schwarzschild BH is the energetically preferred one. 

The aim of this paper is to perform a fully non-linear 
numerical evolution of the scalar solitons, preparing the 
way for a future analysis of the scalar- hairy black holes. 
The philosophy of our analysis is similar to the one of 
Straumman and Zhou for the case of "colored solitons" 
(solitons in Einstein- Yang-Mills theory) The initial 
conditions correspond to unstable scalar solitons in a 
globally regular space-time. Two different sets of ini- 
tial perturbations will be considered: One that leads to 
the formation of a Schwarzschild BH accompanied with 
a small amount of radiated scalar field, and another one 



which corresponds to an "exploding" configuration where 
a global phase transition is triggered through the forma- 
tion of an outward moving domain wall. 

This paper is organized as follows. In Section ^] we 
derive the system of evolution equations and constraints. 
Section ITTTl describes the initial conditions corresponding 
to a static soliton with a small perturbation. We discuss 
the slicing conditions we use for our simulations in Sec- 
tion IIVI Section describes our numerical techniques. 
In Section IVIl we show the numerical results for the two 
types of perturbations mentioned above. We conclude in 
Section IVlII In the Appendix we discuss the hyperbolic- 
ity properties of our system of evolution equation. 



II. FIELD EQUATIONS 

We will consider a model of a scalar field minimally 
coupled to gravity and with a non-trivial self-interaction 
potential. The model is described by the Lagrangian (we 
will use units such that G = c = 1): 



C = 



lOTT 2 



(2.1) 



We choose the following asymmetric scalar-field poten- 
tial leading to the desired asymptotically flat solutions: 



F(0) = -(0-a)2 
4(771 + 772) , 



(2.2) 
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with (T, ?7i, 772 and a constant parameters. For this class 
of potentials one can easily show that, if 771 > 2?72 > 0, 
(j) — a corresponds to a local minimum, = a + 771 to the 
global minimum and ^ = a-l-772 to a local maximum. The 
key feature of this potential for the asymptotically flat 
and static solutions to exist is that the local minimum 
at = a is also a zero of V{(j)) The factor a in 
front of the potential fixes the scale, so one can always 
take (T 1 and just re-scale everything for a different a 
afterward. For the simulations discussed here wc will take 
the following values for the parameters: a — 1, rji = 0.5, 
772 = 0.1 and a = (see Figure P). 
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FIG. 1: Scalar potential V{<l)) corresponding to Eq. 12.21 with 
parameters a = 1, rji — 0.5, 772 =0.1 and a = 0. 



We shall focus on the dynamics of a spherically sym- 
metric space-time described by the metric 



Adr^ + Br^ {dO^ + sin^ dd^'^) . (2.13) 



The spherical symmetry implies that all dynamical func- 
tions depend only on r and t. 

To write down our evolution equations we start by 
defining the quantities 



(2.14) 
(2.15) 
(2.16) 



We will work with the mixed components of the ex- 
trinsic curvature Ka '■= , Kb '■= Kg = K'^, and 
with the matter variables J a Jr, Ma '■— and 
Mb ■— Mg = M^. We also introduce the extra variables 



Da 


:= dr\nA, 


Db 


:= drlnB , 




:= dr In a . 



The field equations following from the Lagrangian (|2.1() 
are the Einstein's field equations and the the Klein- 
Gordon (KG) equation: 



□ c 



dvi4>) 



The stress-energy tensor for the scalar field is 



iv„0V"0 + y(0) 



(2.3) 



(2.4) 



In order to perform a numerical analysis of the prob- 
lem at hand we shall use a 3-1-1 approach based on the 
standard ADM equations Moreover, we shall as- 

sume that the shift vanishes. The evolution equations for 
the 3-metric (7ij) and the extrinsic curvature {Kij) are 



dtK.j = -V,Vja + a{R,j+KK,j 



(2.5) 



2KaK\ ~ SnA'U 



, (2.6) 

and the Hamiltonian and momentum constraints are 

n := R + K^ - K,jK'^ -IGnp^O , (2.7) 
M,; := Vi {K\ - K6\) - Stt J, ^ , (2.8) 

with a the lapse function, 'Di and Rij the covari- 
ant derivative and Ricci tensor associated with 7^, 
R := tvRij, K := tvKij, and where the matter sources 
are defined in terms of the stress-energy tensor as 



T — —77 

Sij — Tij , 



13 = + - (p - 5) 



(2.9) 
(2.10) 
(2.11) 

(2.12) 



with the unit normal to the spatial hypersurfaces. 



D := Da- 2Db , (2.17) 
K := trK = Ka + 2Kb , (2.18) 

and use them instead of Da and Ka-^ 

The evolution equations for {A, B, D, Db, K, Kb} can 
be obtained directly from the ADM equations. However, 
in order to have a hyperbolic evolution system (see Ap- 
pendix) we will remove the terms proportional to OtDb 
and drKs from the ADM evolution equations for K and 
D, respectively, using the constraints. 

Spherical coordinates can be problematic at the ori- 
gin. In order to regularize the coordinate singularity we 
use the regularization scheme described in Ref. [l^] and 
introduce the auxiliary variable 



A 



1 - 



(2.19) 



which will be promoted to an independent dynamical 
quantity and evolved explicitly in time. 

Since we have used the momentum constraint to mod- 
ify the evolution equation for D, for regularizing the 
equations we need to replace this variable with 



ABX 
~A 



The final set of dynamical variables is then 

{A,B,\U,Db,K,Kb\ , 

and their (regularized) evolution equations are 

dtA = 2aA {2Kb - K) , 
dtB = -2aBKB , 
2aA 



B 



drKB + ^ttJa 



^(K-SKb) 



(2.20) 



(2.21) 



(2.22) 
(2.23) 



(2.24) 
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dtU 



-2dr{aK) + 4aDB{K - 3Kb) 



+ 8a 
dtDB = - 
dtK = a 
D. 



DaKB + -^{SKb -K)- in. J a 



2dr{aKB) 

drDa 
A 



(2.25) 
(2.26) 



+ 21'^ 



- AKKb + &Kl + K"^ 
Ar 



dtKB = a 



A J A 

Sir {Ma + 2Mb -2p) 

drDB D^Db , Db 



(2.27) 



2A 2A 

KKb - 8ttMb 
U 2\B 



AA 



U 



a 
Ar 



A 



(2.28) 



These equations can be easily shown to form a strongly 
hyperbolic system (see Appendix). 

In terms of the new variables, the Hamiltonian and 
momentum constraints take the form 

= W := drDB -'^{u -\-Db + 

- AKB{2K-iKB)-DB(^ + ^+^-^) 



+ %ttAp , 
= M:= drKB -{K- 3Kb) 

and the matter variables become 
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Ma = ^ + 
Mb ^ y ■ 

Finally, the KG equation can be written as 

dt'\> = 



an 



BJA 



dti = dr 



an 
bVa 



dtU = -aBy/A d^V + ^dr ( 



where we have introduced the variables 

B^ 



n 



dv 



(2.29) 
(2.30) 

(2.31) 

(2.32) 

(2.33) 
(2.34) 

(2.35) 
(2.36) 
(2.37) 

(2.38) 



For numerical purposes, the evolution equation for H 
above is further transformed into the equivalent form: 



StH = -aB\/A dsV + 3 



aBr^^ 



A 



(2.39) 



The last term on the right hand side of this equation 
includes a first derivative with respect to r^. The reason 
behind this transformation is related to the regularization 
near the origin of the l/r^ factor in the original equation. 
Our final system of evolution equations is then ()2.22|l - 

(E2H1), ESS, and (E23). 



III. INITIAL CONDITIONS 

The initial conditions used to study the evolution 
are the static soliton solutions computed in Ref. jj], 
plus/minus a small Gaussian perturbation in 11. Strictly 
speaking, pure static initial conditions would not result in 
non-trivial evolution. However, since the configurations 
are unstable, the truncation errors are in fact enough to 
trigger the evolution. Nonetheless, since the numerical 
errors become smaller when resolution is increased, it is 
better to include a finite perturbation as "detonator" . 



A. Static soliton 

As mentioned above, we shall consider small pertur- 
bations to the the static soliton found in Ref. IJ. The 
static soliton is obtained by first taking 



B = 1 



Db = Kb=K = 0, 

n = . 



and then solving the set of equations 
drA = A 



1-/4 / 

+ SirAr i^ + V 

r \2A 



dr^ = 


< 


D 


dj-a ~ 


a 




^2^ 


dr(f> — 







Da 
2 



Ad^V , 



(3.1) 
(3.2) 



(3.3) 

(3.4) 

(3.5) 
(3.6) 



for A, a and (j). The first of these equations follows from 
the Hamiltonian constraint (|2.29(l . the second from the 
KG equation (|2.37|) and the last is just the definition of 
£. The third equation is a combination of the evolution 
equation for Kb, Eq. (|2.28|l . and the Hamiltonian con- 
straint. It The equations above are solved subject to the 
boundary conditions 

a{r = 0) = ao , A(r = 0) = 1 , (3.7) 
</.(r = O) = 0o, e(r = 0) = 0. (3.8) 
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FIG. 2; Scalar field (j) for the static configuration (the radial 
coordinate is shown on a logarithmic scale). 
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FIG. 3: Same as Figure|5|for the lapse function a. 



The value (jjo is used as a shooting parameter such that 
the scalar-field settles asymptotically to the local mini- 
mum 1^ = of the potential, thus guaranteeing asymp- 
totic flatness. The value ao, on the other hand, is com- 
pletely arbitrary. The reason for this is that a appears 
only in two of the equations (in one of them through the 
combination — dj-a/a) and both these equations are 
in fact linear in a. In practice one takes ao = 1, solves 
the system of equations, and later re-scales the lapse to 
make sure that its asymptotic value is 1. The numerical 
code we use to solve for the initial data is fourth order 
accurate, in contrast with our evolution code which is 
only second order accurate (see Sect. lVIjl . 

For the particular values of the parameters in the po- 
tential used here, the static soliton solution is shown in 
Figures [21 to 01 In all the plots the numerical grid ex- 
tends to r = 100 and the radial coordinate is shown on a 
logarithmic scale. The value we find for 0o is 



0.40594250807 ±2 x 10 



-11 



(3.9) 



The uncertainty in the eleventh significative figure is es- 
timated by taking the difference from the results of our 
two finest grids with Ar = 0.025 and Ar = 0.0125. It 
is also important to mention that, having fixed the pa- 
rameters of the potential, this static solution is unique, 
as opposed to the case of the corresponding hairy black 
holes where the family of static solutions is parametrized 
by the horizon radius 0]. 



B. Perturbed soliton 

In order to perturb the static soliton we add a small 
Gaussian to the time derivative of the scalar field 



n 



2 / 2 

-r I s 



e << 1 , 



(3.10) 



but still take B = 1, Db = K ^ (notice that it is now 
inconsistent to ask for Kb = as well). We also keep 




log, 0(1 +r) 

FIG. 4: Same as Figure|5|for the metric function A. 



the conditions 9(11 = dtK = 0. The new initial-data set 
of equations to be solved is 



drA = A 



1~A 



SArKl 



SirAr ( + V 



(3.11) 
Ad^V , (3.12) 



dr 



a 



Da 



aA 



6Kl + 87t\^-V 



TP 
A 



drKB = 4^ ^ - 3 — 



(3.13) 
(3.14) 
(3.15) 



where the last equation arises from the fact that the mo- 
mentum constraint is no longer trivial. Notice also that 
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we now have a second order equation for the lapse a. This 
is because we no longer have the condition diKs = 0, but 
rather the condition dtK — 0. In the static case both 
conditions are equivalent, so one can solve the simpler 
first order equation corresponding to dtKs = 0, but in 
the perturbed case this is no longer true. 

The data obtained in this way will be consistent with 
the constraints but will not be static, and will reduce 
to the static soliton solution when e = 0. For the cases 
discussed here, we have taken the Gaussian parameters 
to be e = ±0.002 and s = 10.0. We do not show plots 
for the initial values for (j), a and A, as they look very 
similar to those of the static soliton. 



C. Mass 

To find the mass of the configurations we notice 
that, since the scalar field decays very rapidly with r, 
in the asymptotic region we will have essentially the 
Schwarzschild metric. If we re-parametrize the radial 
metric A as 

(3.16) 

then the mass of the configuration can be obtained from 
M = lim m(r) = lim - f 1 - ^ ) . (3.17) 

r — >oo 1 — ^oo 2 \ A J 

In practice, there is no need to go very far, since in the 
region where the scalar field is zero this expression con- 
verges very rapidly to a limiting value. The mass function 
m{r) for the static configuration is shown in Figure|Sl Its 
limiting value for large r turns out to be 

M = 3.82719754567±2x 10"" , (3.18) 

where again the error bar is estimated by comparing re- 
sults from runs with Ar = 0.025 and Ar = 0.0125. 

For the static configuration one can show, using the 
Hamiltonian constraint, that the mass given by ()3.17fl 
can also be computed as an integral of the energy density 
p associated with the scalar field: 

Af = 47r / pr^dr . (3.19) 
Jo 

However, this is only true for static initial data, as ex- 
trinsic curvature terms will enter this expression for non- 
static configurations. Using the integral expression we 
find the following value for the mass of the static soliton 
(we have only integrated up to the boundary of the grid, 
but since p goes to zero rapidly, the value of the integral 
stops changing to the last decimal figure much before we 
reach the boundary) 

M = 3.827197 ± 2 x lO"'' . (3.20) 



CM - 
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FIG. 5: Same as Figure|5|for the mass function m(r). 

The error in this value is considerably larger, which is not 
surprising since the numerical integral was done using 
the second order trapezium rule. Still, it agrees with the 
result above to 7 significant figures. 

For the perturbed configurations, the value of M com- 
puted with (|3.17l) turns out to be extremely close to the 
value for the static soliton, differing by about 0.1%: 

Mperturbcd = 3.8308619778 ± 8 x 10"^° . (3.21) 

Notice that the mass is the same for both perturbations, 
which is not surprising as p depends only on 11^. 



IV. SLICING CONDITIONS 

We have assumed so far a vanishing shift. We now 
need to specify a slicing condition, that is, a way to find 
the value of the lapse function a during the evolution. A 
common choice is the so-called polar slicing, which in the 
case of vanishing shift is equivalent to the choice of areal 
(or radial) coordinates throughout the evolution. This 
consists of imposing B{t,r) = 1, which in turn implies 
KB{t,r) — 0, and hence K — Ka- From the evolution 
equation for Kb one then finds an ordinary differential 
equation in r for the lapse that can be solved at each 
time step. The main drawback of this approach is that 
the radial-polar slicing gauge does not penetrate BH hori- 
zons, since inside a horizon it is impossible to keep the 
areas of spheres fixed without a non-trivial shift vector. 

Since we will be interested in looking for BH horizons 
during the evolution, we have decided to use horizon 
penetrating coordinates instead of the above gauge. In 
fact, we will use two different types of slicing conditions, 
namely, harmonic slicing and maximal slicing (see be- 
low), each adapted to the physical situation that results 
from a given type of perturbation. The adequate choice, 
of course, is only known a posteriori, so in practice we 
have performed short trial runs in each case with har- 
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monic slicing in order to gain some insight about the 
dynamic behavior of the system. 

Harmonic shcing is a well known condition that re- 
lates the lapse function to the spatial volume elements: 
a = /(a;') 7^/^, with 7 the determinant of the spatial 
metric and / (a:*) an arbitrary time-independent function. 
This slicing condition is in fact equivalent to the require- 
ment that the time coordinate t satisfies the wave equa- 
tion Dt = 0, i.e. i is a harmonic function. The harmonic 
condition can also be written as an evolution equation 
for the lapse in the form 



(4.1) 



which is a particular case of the more general Bona-Masso 
family of slicing conditions pj| . 

It can be shown that harmonic slicing avoids so-called 
"focusing singularities" 0, 0| (those for which the spa- 
tial volume elements vanish at a bounded rate): the lapse 
collapses to zero at the same rate as the volume elements. 
The singularity avoidance of harmonic slicing, however, 
is only marginal in the sense that the slices come arbi- 
trarily close to the singularity after a finite time. For this 
reason harmonic slicing is usually not used for black hole 
evolutions where one wishes to avoid the singularity. 

For black hole evolutions, choices different than har- 
monic are usually better. In particular, maximal slicing, 
defined hy K = dtK — Q, has been a standard workhorse 
for evolving black holes because of its strong singular- 
ity avoidance properties (see e.g. 0|). Maximal slicing 
leads to the following elliptic equation for the lapse 



V^a = a [KijK'^ + An {p + S)] 



(4.2) 



In our case, we have chosen initial data that has precisely 
the property that K = dfK = (even in the perturbed 
case). Moreover, as we will show below, one particular 
perturbation of the scalar soliton leads to a BH forma- 
tion, so maximal slicing is a natural choice. 

In spherical symmetry, the maximal slicing condition 
reduces to the following second order ordinary differential 
equation for the lapse function a 



+4^ ip + SA + 2SB)]. 



(4.3) 



This equation is solved by imposing the following bound- 
ary conditions 



dra — , 

r=Q 



dra =(l-a)/r, (4.4) 

r=ri, 



with ri, the position of the outer boundary. The first 
condition is just the regularity condition at the origin, 
and the second is a Robin outer boundary condition 
that guarantees that the lapse behaves asymptotically 
as 1 + k/r (with k constant). 

Now, although maximal slicing is well adapted for the 
case when the scalar field collapses to a BH, for the 



"exploding" configuration (i.e. an outward moving do- 
main wall) , the harmonic slicing condition is much better 
suited l|4.HI . This is because, as shown in Sec. IVII below, 
in that case the spacetime inside the outward moving wall 
behaves in a manner similar to an anti-de-Sitter space- 
time, as the scalar field moves toward the true minimum 
of the potential. This behavior will produce a big-crunch 
type singularity in a finite proper time, to which maximal 
slicing responds by making the lapse collapse extremely 
rapidly everywhere, thus completely freezing the evolu- 
tion [isIITgI I. Harmonic slicing, on the other hand, makes 
the lapse collapse much slower, allowing a much longer 
evolution. As mentioned above, such slicing allows the 
hypersurfaces to move arbitrarily close to the singularity. 



V. NUMERICAL METHODOLOGY 

For the time integration in our code we use an iterative 
Crank-Nicholson scheme with 3 iterations (see e.g. flTj). 
Derivatives are represented by second order centered fi- 
nite differences on the radial grid. The numerical evolu- 
tion is therefore expected to be second order accurate in 
both Ar and At. We also typically take At = Ar/4 in 
order to be sure that we satisfy the Courant-Friedrichs- 
Levy stability condition. 



A. Boundary conditions 



From the evolution equations (|2.22|) . (|2.23|l and 

it is clear that the metric functions A and B and the 
scalar field (f) evolve only through source terms and can 
be updated point-wise all the way to the boundary of the 
numerical grid. The evolution equations for other vari- 
ables, however, have spatial derivatives in the right hand 
side so we require a boundary condition. For these vari- 
ables we use an outgoing wave (Sommerfeld) boundary 
condition. That is, we assume that near the boundary 
all dynamical variables behave as spherical waves: 



u{r — vt) 



(5.1) 



where v is the wave speed. In practice, we do not use the 
boundary condition H5.1|l as above, but rather we use it 
in differential form 



dth + vdrh + V— — 



(5.2) 



We also assume that the boundary is sufficiently far from 
the origin for the wave speed to be essentially unit (a and 
A are close to 1). So we actually use 



dth — --drh 

r 



(5.3) 



We use finite differences for the above equation (consis- 
tent to second order in both space and time) and apply it 
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to {A, [/, Db, K, Kb,£,, n}. When using harmonic shcing, 
we apply the condition also to and evolve a point- 
wise all the way to the boundary. 

It is important to mention that the boundary con- 
ditions just described are in fact not compatible with 
the constraints. We therefore expect a small amount of 
constraint violation to be introduced at the boundaries, 
which will later propagate inward. We have checked that 
this has essentially no impact on the results presented 
here by moving the outer boundary to different locations 
and comparing the results. Constraint preserving bound- 
ary conditions are of course possible to implement, and 
are certainly desirable, but they are not necessary for the 
results discussed here. 

B. Apparent Horizon 

Throughout the evolution we look for apparent hori- 
zons as an indicator of the formation of a black hole. Of 
course, an apparent horizon will only coincide with an 
event horizon if the spacetime reaches a stationary state 
at late times. But if an apparent horizon is present, an 
event horizon is guaranteed to exist outside of it as long 
as the cosmic censor conjecture and the null energy condi- 
tion hold 18, 19]. Apparent horizons have the advantage 
that their definition is local in the sense that they can 
be located within a given spatial hypersurface. Event 
horizons, on the other hand, are defined globally and can 
therefore only be located once the whole evolution of the 
spacetime is known (or at least, once the evolution is 
known up to the point where the spacetime is essentially 
stationary). In our case, the null energy condition is in- 
dependent of V{(j)), and will therefore hold since for a null 
vector £^ we will have T^^^^^"" = (^''V^c/))^ > 0. More- 
over, as it turns out, when an apparent horizon forms in 
our simulations the scalar field is in the region where the 
potential is positive, so an event horizon is guaranteed to 
exist out of it. 

An apparent horizon (AH) is defined as the outermost 
marginally trapped surface 18], that is, a closed two- 
dimensional surface for which the expansion of the out- 
going null geodesies is zero. In terms of 3-1-1 quantities, 
the expansion of a congruence of null rays moving in 
the outward normal direction to a given surface takes the 
form 

H = Vis' + K,js's' - K , (5.4) 

with the unit outward pointing normal vector to the 
surface. An AH is then the outer-most closed surface 
such that H = everywhere on the surface. 

In spherical symmetry, this expression can be shown 
to reduce to the simple form 

In the code we evaluate this expression over the whole nu- 
merical grid looking for places where it becomes negative 



(notice that for Minkowski we simply have H = 2/r > 0). 
If H ever changes sign, the outermost place where this 
happens is identified as the AH. 

Once an apparent horizon has been located at r = rAH, 
its area will be given by: 

Aah ^ ^TTr\fjBAH (5.6) 

where Bah is the metric function B evaluated at 
r = rAH- We can also associate a mass to this AH 
through the formula 



Mah = v/AWOM = '^"Y ■ (5.7) 

VI. RESULTS FROM NUMERICAL 
SIMULATIONS 

As already mentioned, the perturbations correspond- 
ing to e = 0.002 and e — —0.002 result in very different 
evolutions, so we will consider each case in turn. 

A. Collapse to a black hole 

We will first consider the perturbed configuration cor- 
responding to e = —0.002. In this case we will use max- 
imal slicing to determine the lapse. We have done nu- 
merical simulations using four different grid resolutions, 
Ar = (0.1, 0.05, 0.025, 0.0125). In aU cases the numerical 
grid extended to r — 100. 

Figures El and Ul show the evolution of the metric func- 
tions A and B. In all the plots solid lines correspond 
to initial and final configurations, and dotted lines to in- 
termediate stages (the separation in time between these 
lines is At = 25). At first, there are small oscillations 
around the initial value, but later on the radial metric 
starts to grow in a form characteristic of the slice stretch- 
ing associated with BH spacetimes. Figure |H1 shows the 
behavior of the lapse function. Again, we see the char- 
acteristic collapse of the lapse indicative of the approach 
to a singularity. 

The corresponding evolution of the scalar field can be 
seen in Figure [HI Notice how the scalar field moves to- 
ward the local minimum at = everywhere. At late 
times, the evolution of the inner regions is frozen due to 
the collapse of the lapse there. By that time, however, 
the scalar field has values below 0.1 everywhere, which 
implies that we are in the region where the potential is 
positive. 

The slice stretching and collapse of the lapse are indica- 
tive of the presence of a BH, but they are not enough to 
prove that one is present, so we have looked for the ap- 
pearance of an AH. Figure ITUl shows the mass associated 
with the AH's found during the simulation. From the 
figure one can see that an AH first appears at t ~ 115. 
The mass associated with this horizon grows slightly at 
first as scalar field falls into the BH, but later settles. 
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FIG. 6: Evolution of the metric function A for the pertur- 
bation with e = —0.002. Notice how at late times the metric 
function grows indicating slice stretching. The circles show 
the location of the apparent horizon. 
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FIG. 7: Evolution of the metric function B for the pertur- 
bation with e = —0.002. The circles show the location of the 
apparent horizon. 



For consistency, this mass should always remain below 
the initial ADM mass of the spacetime. From the figure 
we see that for low resolution the numerical error causes 
the horizon mass to eventually become larger than the 
ADM mass. At higher resolutions, however, the horizon 
mass remains below the ADM mass throughout the en- 
tire evolution (see inset of Fig. I1U|) . The small difference 
between the horizon mass and the ADM mass indicates 
that a small amount of scalar field has been radiated 
away. 

A crucial test of the validity of our code is the behav- 
ior of the constraints. Analytically, the constraints must 
be identically zero, but numerical errors mean that for 
our simulations the constraints have in fact finite values. 
Still, those values should converge to zero as resolution 
is increased. Figure ITTI shows the logarithm of the root 




log,o(1+r) 



FIG. 8: Evolution of the lapse function a for the perturbation 
with e = —0.002. Notice the collapse of the lapse at late times, 
indicative of the approach to a singularity. The circles show 
the location of the apparent horizon. 
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FIG. 9: Evolution of the scalar field for the perturbation 
with e = —0.002. At late times the scalar field has values 
below 0.1 everywhere, which implies that we are in the region 
where the potential is positive. 



mean square of the Hamiltonian constraint over the nu- 
merical grid as a function of time for the four resolutions 
used in our simulations Ar = (0.1,0.05,0.025,0.0125). 
As expected, the Hamiltonian constraint is converging to 
zero to second order in Ar. 



B. Explosion 

Next we consider the perturbed configuration corre- 
sponding to e = +0.002. In this case we will use harmonic 
slicing to determine the lapse. We will use the same grid 
resolutions as before, but will extend the numerical grid 
further out to r = 500. 

Figures El and E| show the evolution of the metric 
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FIG. 10: Apparent horizon mass for the perturbation with 
e = -0.002. An AH first appears at t ^ 115. The dashed 
line indicates the initial ADM mass of the spacetime, and the 
solid and dash-dotted lines the horizon mass for the highest 
and lowest resolutions respectively. The asterisks mark the 
first appearance of the apparent horizon. 
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FIG. 11: Logarithm of the root mean square of the Hamil- 
tonian constraint as a function of time for the perturbation 
with e = —0.002, for four difi'erent resolutions. 



functions A and B for this case (lines are now separated 
in time by At = 75). It is evident that the dynamical 
evolution is completely different to the case described in 
the previous section. In the first place, there is no in- 
dication of the slice stretching effect. Moreover, in the 
evolution of the radial metric A it is clear that there is 
a wall moving outward. The wall moves essentially at a 
uniform speed, even if this is not evident in the log plot 
(this speed is approximately 1 in our units, which coin- 
cides with the speed of light in the outer regions). Inside 
this wall the radial metric is collapsing to zero. The an- 
gular metric is also collapsing to zero in this region, but 
not as rapidly. 

The evolution of the lapse function is shown in Fig. PPil 
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FIG. 12: Evolution of the metric function A for the pertur- 
bation with e = +0.002. 
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FIG. 13: Same as Figure [T^ for the metric function B. 

Again, the presence of an outward moving wall is evident. 
Inside the wall the lapse is collapsing to zero, indicating 
that we are approaching a singularity. However, in this 
case no apparent horizon was located for the duration of 
the run. 

The evolution of the scalar field can be seen in Fig. El 
In contrast to the results of the previous section, in this 
case the scalar field is moving toward the true minimum 
of the potential at = 0.5. Since this minimum corre- 
sponds to a negative value of the potential, the interior 
of the wall resembles an anti-de-Sitter spacetime, except 
for the fact that the scalar field is not uniform. Still, 
one would expect the formation of a big-crunch type sin- 
gularity in this region in a finite proper time Il6j |. 
However, because of the singularity avoiding properties of 
harmonic slicing, this singularity would only be reached 
after an infinite coordinate time. 

Finally, Figure El shows the evolution of the logarithm 
of the root mean square of the Hamiltonian constraint for 
the same four resolutions. Second order convergence is 
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FIG. 14: Same as Figure IT^ for the lapse function a. 
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FIG. 15: Evolution of the scalar field 4> for the perturbation 
with e = +0.002. The scalar field is now moving toward the 
true minimum of the potential a.t cj> = 0.5 everywhere. 

again evident. 

The simulation seems to indicate that the perturba- 
tion has triggered a global phase transition from the false 
vacuum in the exterior to the true vacuum in the inte- 
rior. This phase transition propagates through a domain 
wall that always moves outward, which implies that the 
spacetime will never reach a stationary state. 

VII. DISCUSSION 
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FIG. 16: Logarithm of the root mean square of the hamil- 
tonian constraint as a function of time for the perturbation 
with e = -1-0.002, for four different resolutions. 



sponse to two different types of perturbations differing in 
the sign of the perturbing term. 

In one case the configuration is found to undergo grav- 
itational collapse, as the scalar field moves toward the 
local minimum of the potential everywhere. A small 
amount of scalar field is radiated away, and the space- 
time finally settles to a stationary Schwarzschild BH. 

For the other type of perturbation the scalar field "ex- 
plodes" , triggering a global phase transition that propa- 
gates through an outward moving domain wall. This wall 
separates an inner region where the scalar field moves 
toward the true minimum of the potential (the true vac- 
uum), and an outer region with the scalar field in the 
local minimum. As the true minimum has V < Q, the in- 
ner region behaves in a manner similar to anti-de-Sitter 
spacetime, and should produce a big-crunch type singu- 
larity in a finite proper time. The spacetime in this case 
does not reach a stationary state, as the domain wall 
always keeps moving outwards. A question arises as to 
whether the singularity that is forming in this case is 
naked, since there is no evidence for an apparent hori- 
zon [i^ [13 • We believe that this is unlikely, as nothing 
special seems to happen during the evolution. More likely 
an event horizon does exist, but this is no BH in the stan- 
dard sense since, first, no trapped surfaces form inside it, 
and second, a stationary state is never reached and the 
BH eventually swallows the whole spacetime. 



We have considered the dynamical evolution of static 
self- gravitating scalar solitons. These soliton configura- 
tions arise for self-interaction potentials ¥{(()) that have 
a local minimum for which V = Q, and a global minimum 
such that V < Q. This allows for static, asymptotically 
flat solutions for which the scalar field interpolates be- 
tween the two minima. The static configurations were 
shown to be unstable jj, and we have studied their re- 
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APPENDIX 

Recently, the need for having a hyperbolic system of 
evolution equations has been stressed by several authors 
as a necessary condition for the well posedness of the 
Cauchy initial data problem (and also for the stability 
of the numerical evolution). In order to comply with 
this requirement we have re-written the ADM equations 
together with the harmonic slicing condition as a first 
order system of the form 



dtw + drW = S , 
with w a first order variables vector 



w := 



(A.1) 



(A.2) 



IM'" the matrix 



/ 

a/2A 







2a 







\ 

a 
a/A 
2a 



(A.3) 



V -2aA/B 0/ 



and S source terms that include no derivatives (the vari- 
ables {a, A, B} evolve only through source terms and 
need not be considered for the hyperbolicity analysis). 
The eigenvalues of the matrix IM'' turn out to be 



a 



(A.4) 



A3,4 = ± 



a 



VA 



•^5,1 



0, 



with corresponding eigenvectors 



6*1,2 = 


B 


6*3,4 = 






6*5 = 


(0,0,0,0,1,0 


6*6 = 


fo, 0,0, 0,0,1 




(A.5) 
(A.6) 



(A.7) 

(A.8) 
(A.9) 
(A.10) 



Clearly, all eigenvalues are real. Also, the eigenvector 
matrix has determinant equal to 5^/2^3 > 0, showing 

the linear independence of the eigenvectors. The system 
of equations is therefore strongly hyperbolic. 

The above analysis applies to the case of harmonic 
slicing. For maximal slicing the argument is similar, 
except that in that case K is no longer a dynamical 
variable (K = 0), and the lapse is obtained from an 
elliptic equation. The reduced system of equations for 
{Db, Kb,U, X} is still strongly hyperbolic, and the full 
system now has both an elliptic and a hyperbolic sector. 
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